Insights form the National Inpatient Sample
knitr::opts_chunk$set(comment = NA)
library(janitor)
library(broom)
library(haven)
library(naniar)
library(mice)
library(kableExtra)
library(tableone)
library(skimr)
library(Matching)
library(cobalt)
library(survival)
library(twang)
library(survey)
library(conflicted)
library(knitr)
library(tidyverse)
conflicts_prefer(dplyr::select(), dplyr::filter())
theme_set(theme_bw()) Emphysema represents distortion in the lung parenchyma that occurs in patients with chronic obstructive pulmonary disease (COPD) and is pathologically defined as permanent enlargement of the alveoli along with destruction of the alveolar walls. Emphysema can lead to lung hyperinflation and air trapping, which contributes to refractory dyspnea. In carefully selected patients with emphysema, lung volume reduction could improve breathlessness, lung function, quality of life and exercise capacity by decreasing lung hyperinflation 1.
Two mechanisms exist for lung volume reduction, the first, lung volume reduction surgery (LVRS), entails wedge resection of the emphysematous tissue 2. Surgical morbidity can be significant, and certain comorbidities can preclude treatment in many patients. The second is bronchoscopic lung volume reduction (BLVR) with endobronchial valves (EBV). It is a procedure where one-way valves are placed into airway segments in the most diseased lung lobe, allowing air to escape the target area but not allowing it to go back in 3. This mechanism allows the eventual deflation of the target lung lobe and the offloading diaphragmatic stretch.
Two Large clinical trials 4,5 resulted in the FDA approval of endobronchial valves for bronchoscopic lung volume reduction. Overall, respiratory complications occur in 31–35% of patients 6. The most common adverse event associated with EBV placement is pneumothorax (approximately 25% of the cases) 7,8, particularly in the immediate post operative period. Other complications include respiratory failure, acute exacerbation of COPD (AECOPD), pneumonia and hemoptysis. The frequency of these complications is significantly lower than pneumothorax in the published clinical trials.
Since FDA approval in 2018, several centers in the United States have been offering bronchoscopic lung volume reduction for patient with severe COPD. Some centers have published individual reports of complication rates 9,10. To date, there has been no study evaluating complication rates post bronchoscopic lung volume reduction on a population level. While pneumothorax is an expected complication of BLVR, acute respiratory failure can lead to dramatic instability and if progressive, can result in endobronchial valve removal.
To evaluate the odds of acute respiratory failure (ARF) and pneumothorax (PTX) in adult patients with a principal diagnosis of emphysema who are undergoing elective bronchoscopic lung volume reduction (BLVR), in comparison to patients with emphysema who are undergoing surgical lung volume reduction (LVRS) in the immediate post operative period
Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural acute respiratory failure in adult patients with emphysema compared to lung volume reduction surgery?
Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural pneumothorax in adult patients with emphysema compared to lung volume reduction surgery?
Data were obtained from the national inpatient sample (NIS); a US administrative database published under the healthcare utilization project (HCUP) by the Agency for Healthcare Research and Quality 11 through a data use agreement.
The NIS is the largest publicly available all-payer inpatient care database in the United States, containing data on more than seven million hospital stays from nearly 1,000 U.S. hospitals including non-federal public hospitals and academic medical centers. The data include discharge-level records on demographics, diagnoses, procedures, length of hospital stay (LOS), and resource utilization. The NIS uses the International Classification of Diseases, Clinical Modification 10 (ICD-10) to code for diagnoses and procedures.
The data were de-identified, not requiring an IRB approval
Adult patients with a primary diagnosis of emphysema admitted electively for lung volume reduction (surgical or bronchoscopic).
The inclusion and exclusion criteria are aimed at ensuring that the patient was solely admitted for lung volume reduction, and that surgery or endobronchial valve placement was not performed for persistent air leak.
Inclusion criteria:
Exclusion criteria
The sample size consists of 518 adult patients with 253 (49%) patients undergoing bronchoscopic lung volume reduction and 265 patients undergoing surgical lung volume reduction (51%).
The study is a retrospective cohort study, two groups of patients with the same baseline disease (emphysema), one group is receiving the intervention of choice (BLVR), and followed during admission for the development of the outcomes after the intervention of choice.
Bronchoscopic lung volume reduction (BLVR) is the intervention of interest. 253 (49%) of the participants out of 518 in my sample have received BLVR and will be compared to surgical lung volume reduction (LVRS), in 265 (51%) of the participants.
The intervention (BLVR) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BH38GZ, 0BH88GZ, 0BH48GZ, 0BH58GZ, 0BH68GZ, 0BH78GZ, 0BH8GZ, 0BH98GZ.
The comparison (LVRS) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BBC0ZZ, 0BBC4ZZ, 0BBD0ZZ, 0BBD4ZZ, 0BBF0ZZ, 0BBF4ZZ, 0BBG0ZZ, 0BBG4ZZ, 0BBH0ZZ, 0BBH4ZZ, 0BBJ0ZZ, 0BBJ4ZZ, 0BBK0ZZ, 0BBK4ZZ, 0BBL0ZZ, 0BBL4ZZ, 0BBM0ZZ, 0BBM4ZZ
The propensity score model used covariates that are related to patient characteristics (age, sex, race, presence of chronic respiratory failure, insurance carrier and classification of the median income for the patient’s zip code) and hospital characteristics (hospital size, location, teaching status, hospital ownership and year of the procedure)
The covariates were measured and recorded on admission and prior to the procedure.
Acute respiratory failure is a new (or increased) need for oxygen supplementation in order to maintain adequate oxygenation.
It will be defined using ICD-10 diagnosis codes J9600, J9601, J9602, J9620, J9621, J9622, J95821, J95822 and is a binary outcome (Yes/No)
Pneumothorax is the abnormal collection of air around the lung leading to lung collapse, and possibly causing respiratory failure or in severe forms, instability and death.
This will be defined using ICD-10 diagnosis codes J930, J9311, J9312, J9383, J939, J95811, J95812 and is a binary outcome (Yes/No)
The final analytical data was prepared using SAS. This included importing the national inpatient sample (NIS) into a statistical program that can handle the large size and number of original observations (roughly 21 million observations). SAS was used to properly create the variables of interest based on ICD-10 codes. A final dataset with the intended variables was then created for the final analyses.
The data was imported in the format of sas7bdat using haven package. zap_labels and clean_names were used to remove column labels and standardize names respectively into R format.
blvr <- read_sas("blvr.sas7bdat", NULL) |> zap_label() |> clean_names() glimpse(blvr)Rows: 540
Columns: 18
$ age <dbl> 70, 69, 73, 57, 67, 55, 44, 55, 65, 43, 81, 69, 63, 7…
$ elective <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ female <dbl> 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0,…
$ pay1 <dbl> 1, 3, 1, 1, 3, 1, 3, 3, 1, 3, 1, 1, 3, 1, 3, 2, 2, 3,…
$ race <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6, 1, 2, 1,…
$ year <dbl> 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018,…
$ zipinc_qrtl <dbl> 2, 2, 4, 2, 2, 4, 2, 4, 4, 4, 3, 4, 4, 3, 3, 3, 1, 1,…
$ hosp_bedsize <dbl> 2, 3, 3, 3, 2, 2, 3, 3, 3, 2, 1, 3, 3, 3, 1, 1, 3, 2,…
$ hosp_locteach <dbl> 3, 3, 3, 3, 3, 3, 1, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3,…
$ hosp_region <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ h_contrl <dbl> 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2,…
$ emphysema <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ ebv <dbl> 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ lvrs <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,…
$ crf <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ ptx <dbl> 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,…
$ arf <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,…
$ volume_reduction <chr> "LVRS", "BLVR", "LVRS", "BLVR", "LVRS", "LVRS", "LVRS…The code below achieved the following steps
ebv to vol_reduction, this was intended to be the numeric 0/1 version of the interventionvol_reduction which is the intervention of choice as a numeric variable (0/1) for further analysisblvr <- blvr |>
filter(age > 17) |> rename(vol_reduction = ebv) |>
mutate(sex = fct_recode(factor(female), "Male" = "0", "Female" = "1"),
insurance = fct_recode(factor(pay1), "Medicare" = "1", "Medicaid" = "2", "Private" = "3", "Private" = "4", "Private" = "6"),
insurance = fct_relevel(insurance, "Medicare", "Medicaid", "Private"),
race = fct_recode(factor(race), "White" = "1", "Black" = "2", "Other" = "3", "Other" = "4", "Other" = "5", "Other" = "6"),
hospital_size = fct_recode(factor(hosp_bedsize), "Small" = "1", "Medium" = "2", "Large" = "3"),
teaching = fct_recode(factor(hosp_locteach), "Rural" = "1", "Urban nonteaching" = "2", "Urban teaching" = "3"),
region = fct_recode(factor(hosp_region), "Northeast" = "1", "Midwest" = "2", "South" = "3", "West" = "4"),
ownership = fct_recode(factor(h_contrl), "Government" = "1", "Private, non-profit" = "2", "Private, for-profit" = "3"),
income = fct_recode(factor(zipinc_qrtl), "0-25th percentile" = "1", "26th to 50th percentile" = "2", "51st to 75th percentile" = "3", "76th to 100th percentile" = "4"),
crf = fct_recode(factor(crf), "Yes" = "1", "No" = "0"),
ptx_f = fct_recode(factor(ptx), "Yes" = "1", "No" = "0"),
arf_f = fct_recode(factor(arf), "Yes" = "1", "No" = "0"),
vol_reduction_f = as_factor(volume_reduction),
vol_reduction_f = fct_relevel(vol_reduction_f, "LVRS"),
year = as_factor(year)) The final analytical tibble included the covariates used in building the propensity score, as well as both numeric/ factor versions of the intervention and outcomes
blvr <- blvr |> select(age, sex, race, crf, income, insurance, hospital_size, teaching, region, ownership, year, vol_reduction, vol_reduction_f, arf, arf_f, ptx, ptx_f)blvr |> miss_var_summary()# A tibble: 17 × 3
variable n_miss pct_miss
<chr> <int> <dbl>
1 race 14 2.70
2 income 9 1.74
3 age 0 0
4 sex 0 0
5 crf 0 0
6 insurance 0 0
7 hospital_size 0 0
8 teaching 0 0
9 region 0 0
10 ownership 0 0
11 year 0 0
12 vol_reduction 0 0
13 vol_reduction_f 0 0
14 arf 0 0
15 arf_f 0 0
16 ptx 0 0
17 ptx_f 0 0 There were fourteen missing observations in the race variable, and nine observations in the income variable. These constituted 4.4% of the total data, and assumed to be missing at random (MAR).
blvr |> tabyl(income) income n percent valid_percent
0-25th percentile 139 0.26833977 0.2730845
26th to 50th percentile 133 0.25675676 0.2612967
51st to 75th percentile 129 0.24903475 0.2534381
76th to 100th percentile 108 0.20849421 0.2121807
<NA> 9 0.01737452 NAblvr |> tabyl(race) race n percent valid_percent
White 421 0.81274131 0.83531746
Black 46 0.08880309 0.09126984
Other 37 0.07142857 0.07341270
<NA> 14 0.02702703 NAGiven the small sample and the desire to keep as many observations as possible, single imputation using the mice package. Fifteen imputed data sets were created and one dataset was later used as the single imputation tibble.
set.seed(123)
blvr_mice <- mice(blvr, m = 15, printFlag = FALSE)blvr <- complete(blvr_mice, 10) |> tibble()Checking the variables with the missing values after imputation showed preserved ratios in each level in comparison to pre-imputation
blvr |> tabyl(income) income n percent
0-25th percentile 142 0.2741313
26th to 50th percentile 136 0.2625483
51st to 75th percentile 131 0.2528958
76th to 100th percentile 109 0.2104247blvr |> tabyl(race) race n percent
White 432 0.83397683
Black 48 0.09266409
Other 38 0.07335907An ID row was added to identify observations if needed
blvr <- tibble::rowid_to_column(blvr, "ID"); blvr <- blvr |> mutate (ID = as.character(ID))write.csv(blvr, "blvr.csv")| Variable | Type | Description | Role |
|---|---|---|---|
age |
Numeric | Age in years | Co-variate |
sex |
Factor (w/2 levels) | Sex of the subject: Male or Female | Covariate |
race |
Factor (w/3 levels) | White, Black, Other | Covariate |
crf |
Factor (w/2 levels) | Presence of chronic respiratory failure: Yes or No | Covariate |
income |
Factor (w/4 levels) | Quartile classification of the estimated median household income of residents in the patient’s ZIP Code: 0-25th percentile, 26-50th percentile, 51-75th percentile and 76-00th percentile | Covariate |
insurance |
Factor (w/3 levels) | Medicare, Medicaid or Private | Covariate |
hospital_size |
Factor (w/3 levels) | Small, Medium, Large | Covariate |
teaching |
Factor (w/3 levels) | Teaching status of the hospital: Rural, Urban non-teaching and Urban teaching | Covariate |
region |
Factor (w/4 levels) | Region of the hospital: Northeast, Midwest, South and West | Covariate |
ownership |
Factor (w/3 levels) | Ownership of the hospital: Government, Private, non-profit and Private for-profit | Covariate |
year |
Character | Year of the procedure | Covariate |
vol_reduction |
numeric | Type of lung volume reduction: bronchoscopic (1) or surgical (0) | Intervention/Comparator |
vol_reduction_f |
Factor (w/2 levels) | Type of lung volume reduction: bronchoscopic BLVR or surgical LVRS |
Intervention/Comparator |
arf |
Numeric | Presence of acute respiratory failure: 0/1 | Outcome 1 |
arf_f |
Factor (w/2 levels) | Presence of acute respiratory failure: Yes or No | Outcome 1 |
ptx |
Numeric | Presence of pneumothorax: 0/1 | Outcome 1 |
ptx_f |
Factor (w/2 levels) | Presence of pneumothorax: Yes or No | Outcome 2 |
blvr.vars <- c("age", "race", "sex", "crf", "income", "insurance",
"hospital_size", "teaching", "region", "ownership", "year")
blvr.reduction <- c("vol_reduction_f")
blvr_table1 <- CreateTableOne(data = blvr,
vars = blvr.vars,
strata = blvr.reduction,
test = FALSE)print(blvr_table1) |> kable() |> kable_material() Stratified by vol_reduction_f
LVRS BLVR
n 265 253
age (mean (SD)) 56.40 (15.07) 67.88 (8.39)
race (%)
White 204 (77.0) 228 (90.1)
Black 31 (11.7) 17 ( 6.7)
Other 30 (11.3) 8 ( 3.2)
sex = Female (%) 76 (28.7) 110 (43.5)
crf = Yes (%) 34 (12.8) 117 (46.2)
income (%)
0-25th percentile 74 (27.9) 68 (26.9)
26th to 50th percentile 68 (25.7) 68 (26.9)
51st to 75th percentile 70 (26.4) 61 (24.1)
76th to 100th percentile 53 (20.0) 56 (22.1)
insurance (%)
Medicare 108 (40.8) 196 (77.5)
Medicaid 51 (19.2) 11 ( 4.3)
Private 106 (40.0) 46 (18.2)
hospital_size (%)
Small 12 ( 4.5) 27 (10.7)
Medium 55 (20.8) 42 (16.6)
Large 198 (74.7) 184 (72.7)
teaching (%)
Rural 7 ( 2.6) 7 ( 2.8)
Urban nonteaching 20 ( 7.5) 28 (11.1)
Urban teaching 238 (89.8) 218 (86.2)
region (%)
Northeast 63 (23.8) 32 (12.6)
Midwest 79 (29.8) 77 (30.4)
South 78 (29.4) 105 (41.5)
West 45 (17.0) 39 (15.4)
ownership (%)
Government 37 (14.0) 31 (12.3)
Private, non-profit 210 (79.2) 214 (84.6)
Private, for-profit 18 ( 6.8) 8 ( 3.2)
year (%)
2018 85 (32.1) 6 ( 2.4)
2019 97 (36.6) 76 (30.0)
2020 83 (31.3) 171 (67.6) | LVRS | BLVR | |
|---|---|---|
| n | 265 | 253 |
| age (mean (SD)) | 56.40 (15.07) | 67.88 (8.39) |
| race (%) | ||
| White | 204 (77.0) | 228 (90.1) |
| Black | 31 (11.7) | 17 ( 6.7) |
| Other | 30 (11.3) | 8 ( 3.2) |
| sex = Female (%) | 76 (28.7) | 110 (43.5) |
| crf = Yes (%) | 34 (12.8) | 117 (46.2) |
| income (%) | ||
| 0-25th percentile | 74 (27.9) | 68 (26.9) |
| 26th to 50th percentile | 68 (25.7) | 68 (26.9) |
| 51st to 75th percentile | 70 (26.4) | 61 (24.1) |
| 76th to 100th percentile | 53 (20.0) | 56 (22.1) |
| insurance (%) | ||
| Medicare | 108 (40.8) | 196 (77.5) |
| Medicaid | 51 (19.2) | 11 ( 4.3) |
| Private | 106 (40.0) | 46 (18.2) |
| hospital_size (%) | ||
| Small | 12 ( 4.5) | 27 (10.7) |
| Medium | 55 (20.8) | 42 (16.6) |
| Large | 198 (74.7) | 184 (72.7) |
| teaching (%) | ||
| Rural | 7 ( 2.6) | 7 ( 2.8) |
| Urban nonteaching | 20 ( 7.5) | 28 (11.1) |
| Urban teaching | 238 (89.8) | 218 (86.2) |
| region (%) | ||
| Northeast | 63 (23.8) | 32 (12.6) |
| Midwest | 79 (29.8) | 77 (30.4) |
| South | 78 (29.4) | 105 (41.5) |
| West | 45 (17.0) | 39 (15.4) |
| ownership (%) | ||
| Government | 37 (14.0) | 31 (12.3) |
| Private, non-profit | 210 (79.2) | 214 (84.6) |
| Private, for-profit | 18 ( 6.8) | 8 ( 3.2) |
| year (%) | ||
| 2018 | 85 (32.1) | 6 ( 2.4) |
| 2019 | 97 (36.6) | 76 (30.0) |
| 2020 | 83 (31.3) | 171 (67.6) |
Based on table 1, patients receiving bronchoscopic lung volume reduction (BLVR) were noted to be older (mean age of 67 years), with a higher percentage of females (43.5%), likely due to a preference of a less invasive approach compared to the surgical group. BLVR patients also tended to have a higher percentage of chronic respiratory failure, likely due to the ease of valve removal should the patient’s baseline respiratory failure worsen after procedure. The patient population was more white (90%), with medicare as the primary insurance provider (77.5%) which is likely due to the age bracket in this group (65 years and older).
As expected both groups had a similar distribution of median income per zip code, and majority of these procedures were done in large, non for profit and urban teaching hospitals. Lung volume reduction surgery number distribution was steady over the years included (roughly around 30%), however bronchoscopic lung volume reduction numbers significantly increased over the years (starting in 2019 and 2020)
arfggplot(blvr, aes(x = vol_reduction_f, fill = arf_f)) +
geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position = "dodge") +
scale_y_continuous(labels = scales::percent) +
scale_fill_viridis_d(option = "turbo", alpha = 0.5) + labs(title = "Respiratory Failure", x = "Volume Reduction", y = "% of Patients", fill = "Resp. Failure")
blvr |> tabyl(arf_f, vol_reduction_f) |> kable() |> kable_material()| arf_f | LVRS | BLVR |
|---|---|---|
| No | 235 | 209 |
| Yes | 30 | 44 |
The visualization above suggested that BLVR patients tended to suffer a higher percentage of acute respiratory failure
A logistic regression model was fit to assess the odds of respiratory failure based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio
unadjusted_arf <- glm(arf ~ vol_reduction, data=blvr, family=binomial)kable(tidy(unadjusted_arf, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_paper()| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.13 | 0.19 | -10.62 | 0.00 | 0.09 | 0.18 |
| vol_reduction | 1.65 | 0.26 | 1.96 | 0.05 | 1.00 | 2.74 |
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.65 odds of acute respiratory failure (95% CI (1,2.74)). This confidence interval contained 1, thus it was difficult to reliably say that BLVR would meaningfully lead to a higher odds of acute respiratory failure
Overall, it seems that the patients receiving BLVR tended to have higher odds of respiratory failure, even though the odds ratio’s interval contained 1, this odds ratio shows an interesting trend
ptxggplot(blvr, aes(x = vol_reduction_f, fill = ptx_f)) +
geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position = "dodge") +
scale_y_continuous(labels = scales::percent) +
scale_fill_viridis_d(option = "turbo", alpha = 0.5) +
labs(title = "Pneumothorax", x = "Volume Reduction", y = "% of Patients", fill = "PTX")
blvr |> tabyl(ptx_f, vol_reduction_f) |> kable() |> kable_material()| ptx_f | LVRS | BLVR |
|---|---|---|
| No | 159 | 200 |
| Yes | 106 | 53 |
The visualization above suggested that BLVR patients tended to suffer a lower percentage of pneumothorax
A logistic regression model was fit to assess the odds of pneumothorax based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio
unadjusted_ptx <- glm(ptx ~ vol_reduction, data=blvr, family=binomial)kable(tidy(unadjusted_ptx, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_paper()| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.67 | 0.13 | -3.23 | 0 | 0.52 | 0.85 |
| vol_reduction | 0.40 | 0.20 | -4.64 | 0 | 0.27 | 0.58 |
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR had a 0.4 odds of pneumothorax (95% CI (0.27,0.58)), which is lower odds than the subject receiving LVRS, which is plausible as LVRS is more invasive.
A logistic regression model was used to estimate the propensity score. Having BLVR/LVRS was modeled as the outcome, with patient and hospital characteristics as independent variables. Patient characteristics used in the propensity score included age, sex, race, crf, income and insurance. Hospital characteristics included hospital_size, teaching status, region, ownership and year of the procedure. The raw and linear propensity score values were saved to the blvr tibble.
psmodel <- glm(vol_reduction ~ age + sex + race + crf + income + insurance +
hospital_size + teaching + region + ownership + year,
family=binomial, data=blvr)
summary(psmodel)
Call:
glm(formula = vol_reduction ~ age + sex + race + crf + income +
insurance + hospital_size + teaching + region + ownership +
year, family = binomial, data = blvr)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.696344 1.503545 -3.124 0.001787 **
age 0.064113 0.015710 4.081 4.48e-05 ***
sexFemale 0.716995 0.258418 2.775 0.005528 **
raceBlack -0.481109 0.451914 -1.065 0.287055
raceOther -1.521225 0.599723 -2.537 0.011195 *
crfYes 1.294080 0.288002 4.493 7.01e-06 ***
income26th to 50th percentile 0.393352 0.350994 1.121 0.262423
income51st to 75th percentile -0.286941 0.342652 -0.837 0.402361
income76th to 100th percentile 0.026792 0.386403 0.069 0.944722
insuranceMedicaid -1.353403 0.515379 -2.626 0.008639 **
insurancePrivate -0.469448 0.318487 -1.474 0.140484
hospital_sizeMedium -2.400828 0.671416 -3.576 0.000349 ***
hospital_sizeLarge -2.077553 0.616456 -3.370 0.000751 ***
teachingUrban nonteaching 0.678464 0.954277 0.711 0.477102
teachingUrban teaching -0.472681 0.823311 -0.574 0.565885
regionMidwest 0.044408 0.388618 0.114 0.909021
regionSouth 0.668451 0.378964 1.764 0.077750 .
regionWest 0.167103 0.437435 0.382 0.702457
ownershipPrivate, non-profit -0.006272 0.357697 -0.018 0.986011
ownershipPrivate, for-profit -1.513048 0.774833 -1.953 0.050850 .
year2019 2.211839 0.517582 4.273 1.93e-05 ***
year2020 3.281191 0.507832 6.461 1.04e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 717.82 on 517 degrees of freedom
Residual deviance: 436.37 on 496 degrees of freedom
AIC: 480.37
Number of Fisher Scoring iterations: 6blvr$ps <- psmodel$fitted
blvr$linps <- psmodel$linear.predictorsggplot(blvr, aes(x = vol_reduction_f, y = ps, color = vol_reduction_f)) +
geom_boxplot() +
geom_jitter(width = 0.1, alpha = 0.3) +
scale_color_viridis_d(option = "plasma", end = 0.5) +
labs(title = "Boxplot & Jitter", subtitle = "Distribution Of the Raw PS", x = "Type of Procedure", y = "Propensity Scores",
caption = "* Line indicates sample median", col = "Procedure")
ggplot(blvr, aes(x = linps, fill = vol_reduction_f)) +
geom_density(alpha = 0.3) +
scale_fill_viridis_d(option = "turbo") +
labs(title = "Density Plot", subtitle = "For the Linear Propensity Score", x = "Linps", y = "Density",
fill = "Volume Reduction")
blvr |> group_by(vol_reduction_f) |> skim_without_charts(ps, linps)| Name | group_by(blvr, vol_reduct… |
| Number of rows | 518 |
| Number of columns | 20 |
| _______________________ | |
| Column type frequency: | |
| numeric | 2 |
| ________________________ | |
| Group variables | vol_reduction_f |
Variable type: numeric
| skim_variable | vol_reduction_f | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 |
|---|---|---|---|---|---|---|---|---|---|---|
| ps | LVRS | 0 | 1 | 0.27 | 0.26 | 0.00 | 0.03 | 0.18 | 0.47 | 0.93 |
| ps | BLVR | 0 | 1 | 0.72 | 0.24 | 0.03 | 0.60 | 0.80 | 0.91 | 1.00 |
| linps | LVRS | 0 | 1 | -1.81 | 2.14 | -7.45 | -3.37 | -1.53 | -0.13 | 2.52 |
| linps | BLVR | 0 | 1 | 1.29 | 1.53 | -3.61 | 0.40 | 1.38 | 2.27 | 5.64 |
The overlap between the propensity scores in intervention groups was assessed numerically and graphically. The boxplot and the density plots show some overlap in the propensity scores between the intervention and the treatment.
The LVRS group had a lower propensity score mean, with problematic values falling below 0.01. The BLVR group also had some problematic high values over 0.99. Troubleshooting will be addressed at a later step
unadjusted_rubin1 <-
abs(100*(mean(blvr$linps[blvr$vol_reduction==1]) -
mean(blvr$linps[blvr$vol_reduction==0])) /
sd(blvr$linps))
unadjusted_rubin2 <- with(blvr, var(linps[vol_reduction == 1]) / var(linps[vol_reduction == 0]))unadjusted_rubin1; unadjusted_rubin2[1] 127.7055[1] 0.506223Rubin’s rule 1 value was higher than 50%, this meant that running an unadjusted model could not be justified without applying the propensity scores to account for selection bias.
Rubin’s rule 2 was on the lower bound of acceptable value as well, as the variance ratio of the linear propensity score in the BLVR (intervention) group to the variance of the linear propensity score in the LVRS (control) group should be close to 1 (and ideally between 0.8 and 1.25)
blvr |> count(ps >= 0.99); blvr |> count(ps < 0.01)# A tibble: 2 × 2
`ps >= 0.99` n
<lgl> <int>
1 FALSE 514
2 TRUE 4# A tibble: 2 × 2
`ps < 0.01` n
<lgl> <int>
1 FALSE 490
2 TRUE 28There were twenty eight observations with propensity score values below 0.01 and four observations with propensity score values above 0.99. The extreme lows seemed to be present in the LVRS group while the extreme highs were present in the BLVR group
mosaic::favstats(ps ~ vol_reduction_f, data = blvr)Registered S3 method overwritten by 'mosaic':
method from
fortify.SpatialPolygonsDataFrame ggplot2 vol_reduction_f min Q1 median Q3 max
1 LVRS 0.0005808267 0.03308973 0.1778799 0.4680281 0.9256507
2 BLVR 0.0263071482 0.59915004 0.7983214 0.9061054 0.9964521
mean sd n missing
1 0.2679640 0.2584781 265 0
2 0.7193262 0.2402497 253 0The following strategies were attempted in order to resolve the above issue:
vif didn’t reveal problematic colinearity between the variables used in the propensity score modelcar::vif(psmodel) GVIF Df GVIF^(1/(2*Df))
age 1.634643 1 1.278532
sex 1.070972 1 1.034878
race 1.204488 2 1.047612
crf 1.105123 1 1.051248
income 1.428943 3 1.061294
insurance 1.667039 2 1.136283
hospital_size 1.403948 2 1.088523
teaching 1.250028 2 1.057377
region 1.470981 3 1.066435
ownership 1.355067 2 1.078922
year 1.164515 2 1.038810Re-constructing the propensity score model, adding one covariate at a time, while checking which variable drives the propensity score values below 0.1 or above 0.99 showed that each of the following covariates race, year, sex, ownership and crf could independently drop the propensity score below 0.01 or drive it up above 0.99. Excluding these covariates form the model would be problematic in itself
Collapsing multicategorical variables was attempted, these variables included year, race, ownership, insurance, hospital_size, teaching and region. It did not result in the improvement of the propensity score even when each covariate was added individually after collapsing.
Filtering obaservations with problematic propensity scores proved to be the most reasonable approach in this situation. This resulted in a sample with 237 observations on the control intervention (LVRS) and 249 observation on the studied intervention (BLVR). The resulting dataset blvr_matching was used for propensity score matching analyses
blvr_matching <- blvr |> filter(ps <= 0.990 & ps >= 0.010)mosaic::favstats(ps ~ vol_reduction_f, data = blvr_matching) vol_reduction_f min Q1 median Q3 max mean
1 LVRS 0.01020130 0.05651206 0.2244966 0.4944178 0.9256507 0.2991748
2 BLVR 0.02630715 0.58669227 0.7927773 0.9014236 0.9892394 0.7149098
sd n missing
1 0.2558878 237 0
2 0.2396084 249 0A density plot was created to describe the overlap of the linear propensity scores between the two intervention groups after balancing the dataset, the overlap is described graphically in the density plot below.
ggplot(blvr_matching, aes(x = linps, fill = vol_reduction_f)) +
geom_density(alpha = 0.3) +
scale_fill_viridis_d(option = "turbo") +
labs(title = "Density Plot", subtitle = "For the Linear Propensity Score", x = "Linps", y = "Density",
fill = "Volume Reduction")
Propensity score matching was conducted using a 1:1 caliper matching on the linear propensity score without replacement from the Matching package. A default caliper value of 0.2 was selected. The type of matching was chosen based on the similar sample size between the intervention and control groups.
X <- blvr_matching$linps
Tr <- as.logical(blvr_matching$vol_reduction)
match <- Match(Tr = Tr, X = X, M = 1,
caliper = 0.2,
estimand = "ATT", replace = FALSE,
ties = FALSE)Warning in Match(Tr = Tr, X = X, M = 1, caliper = 0.2, estimand = "ATT", :
replace==FALSE, but there are more (weighted) treated obs than control obs.
Some treated obs will not be matched. You may want to estimate ATC instead.summary(match)
Estimate... 0
SE......... 0
T-stat..... NaN
p.val...... NA
Original number of observations.............. 486
Original number of treated obs............... 249
Matched number of observations............... 107
Matched number of observations (unweighted). 107
Caliper (SDs)........................................ 0.2
Number of obs dropped by 'exact' or 'caliper' 142 One hundred and seven BLVR patients were matched to one hundred and seven LVRS patients. The degree of covariate imbalance was assessed before and after propensity score matching for each covariate as well as the linear and raw propensity scores
covs <- blvr_matching |>
select(age, sex, race, crf, insurance, income, hospital_size, region, teaching, ownership, year, ps, linps)
b <- bal.tab(match,
treat = blvr_matching$vol_reduction,
covs = covs,
stats = c("m", "v"),
quick = FALSE, un = TRUE, disp.v.ratio = TRUE)
bBalance Measures
Type Diff.Un V.Ratio.Un Diff.Adj V.Ratio.Adj
age Contin. 1.1098 0.3830 0.1295 0.8005
sex_Female Binary 0.1135 . -0.0280 .
race_White Binary 0.0977 . 0.0093 .
race_Black Binary -0.0497 . 0.0093 .
race_Other Binary -0.0480 . -0.0187 .
crf_Yes Binary 0.3104 . 0.0187 .
insurance_Medicare Binary 0.3238 . 0.0841 .
insurance_Medicaid Binary -0.1204 . -0.0093 .
insurance_Private Binary -0.2034 . -0.0748 .
income_0-25th percentile Binary -0.0052 . 0.0000 .
income_26th to 50th percentile Binary -0.0008 . -0.0280 .
income_51st to 75th percentile Binary -0.0333 . 0.0280 .
income_76th to 100th percentile Binary 0.0392 . 0.0000 .
hospital_size_Small Binary 0.0542 . -0.0187 .
hospital_size_Medium Binary -0.0170 . 0.0000 .
hospital_size_Large Binary -0.0372 . 0.0187 .
region_Northeast Binary -0.0951 . 0.0374 .
region_Midwest Binary 0.0181 . 0.0187 .
region_South Binary 0.0934 . 0.0000 .
region_West Binary -0.0164 . -0.0561 .
teaching_Rural Binary -0.0014 . 0.0093 .
teaching_Urban nonteaching Binary 0.0367 . -0.0093 .
teaching_Urban teaching Binary -0.0353 . 0.0000 .
ownership_Government Binary -0.0314 . -0.0467 .
ownership_Private, non-profit Binary 0.0499 . 0.0280 .
ownership_Private, for-profit Binary -0.0185 . 0.0187 .
year_2018 Binary -0.2333 . -0.0467 .
year_2019 Binary -0.0956 . 0.1402 .
year_2020 Binary 0.3289 . -0.0935 .
ps Contin. 1.7351 0.8768 0.0429 1.0478
linps Contin. 1.7602 0.7249 0.0326 1.0488
Sample sizes
Control Treated
All 237 249
Matched 107 107
Unmatched 130 0
Discarded 0 142bal.plot(match,
treat = blvr_matching$vol_reduction,
covs = covs,
var.name = "linps",
which = "both",
sample.names =
c("Unmatched Sample", "Matched Sample")) +
scale_fill_viridis_d(option = "turbo") +
labs(title = "Distributional Balance for linear PS", subtitle = "1:1 Caliper Matching", x = "Linear PS") 
A Love plot of standardized differences before and after matching was obtained to evaluate matching success, showing improvement over the unmatched sample
love.plot(b, threshold = .1, size = 3,
var.order = "unadjusted",
colors = c("hotpink3", "purple4"),
stars = "raw") +
labs(title = "Love Plot of Standardized Differences", subtitle = "1:1 Caliper Matching")
Rubin’s first and second rules were assessed after matching, showing improvement over the unadjusted Rubin’s rules.
matches <- factor(rep(match$index.treated, 2))
blvr.matchedsample <- cbind(matches, blvr_matching[c(match$index.control, match$index.treated),])matched_rubin1 <- with(blvr.matchedsample,
abs(100*(mean(linps[vol_reduction==1])-mean(linps[vol_reduction==0]))/sd(linps)))
matched_rubin2 <- with(blvr.matchedsample, var(linps[vol_reduction==1])/var(linps[vol_reduction==0]))matched_rubin1; matched_rubin2[1] 4.060997[1] 1.048778Rubin’s rule 1 was less than 10%, and Rubin’s rule 2 was around 1 which was ideal.
arfmatched_arf <- clogit(arf ~ vol_reduction + strata(matches), data=blvr.matchedsample)kable(tidy(matched_arf, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_material()| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| vol_reduction | 1.73 | 0.38 | 1.44 | 0.15 | 0.82 | 3.63 |
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.73 odds of acute respiratory failure (95% CI (0.82,3.63)). The subject receiving BLVR tended to have a higher odds ratio of respiratory failue in comparison to patients receiving LVRS, however this confidence interval crosses 1
ptxmatched_ptx <- clogit(ptx ~ vol_reduction + strata(matches), data=blvr.matchedsample)kable(tidy(matched_ptx, exponentiate = TRUE, conf.level = 0.95, conf.int = TRUE), digits = 2) |> kable_material()| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| vol_reduction | 0.68 | 0.3 | -1.3 | 0.19 | 0.38 | 1.22 |
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 0.68 odds of pneumothorax (95% CI (0.38,1.22)). Patients receiving BLVR seemed to have a higher odds ratio of pneumothorax in comparison to patients receiving BLVR, however this confidence interval contained 1, the results were not statistically detectable and a stability analysis instead of a sensitivity analysis was carried out.
Trying different imputed datasets (generated through mice package) resulted in similar outcome estimates (odds ratio)
The matching strategies described below were attempted. None of these strategies provided an improvement in the effective sample size, propensity score distributional balance, standardized differences or Rubin’s rules in comparison to the matching strategy used in this analysis (1:1 matching with a caliper)
| Match | Trtmnt | Ctrl | Distributional Balance | Std Differences | Rubin’s 1 | Rubin’s 2 |
|---|---|---|---|---|---|---|
1:1 Greedy W replacement |
249 | 75 | Better overlap, not ideal | better but dispersed | 8.8 | 1.27 |
Genetic matching |
273 | 78 | Better overlap, not ideal | Dispersed | 8.8 | 1.26 |
1:1 Neartest neighbor matching W Replacement |
249 | 79 | Better overlap, not ideal | better but dispersed (> 0.1) | 8 | 1.21 |
Given the initial difference in sample sizes across the two volume reduction cohorts, propensity score matching has left a selection of the patients unmatched. Weighting by the inverse propensity score (average of the treatment effect on the treated) was next carried out using the initial database blvr (the database that still contained propensity scores below 0.01 and above 0.99)
blvr$wts1 <- ifelse(blvr$vol_reduction==1, 1, blvr$ps/(1-blvr$ps))Then the resulting ATT (average treatment effect on the treated) weights were plotted in the figure below, where the treatment group got the same weight of 1, and the control group had varying weights
ggplot(blvr, aes(x = ps, y = wts1, color = vol_reduction_f)) +
geom_point(size = 2, alpha = 0.4) +
facet_wrap(~ vol_reduction_f) +
labs(x = "Estimated Propensity for Volume Reduction",
y = "ATT weights for blvr",
title = "ATT weighting structure",
subtitle = "In the BLVR dataset",
col = "Volume Reduction")
The balance after weighting was assessed using dx.wts from the twang package. Prior to weighting, the control sample was 237 participants, and after weighting, the effective sample size was 54 participants, representing a sample size that would be required to achieve the same level of precision if that sample were a simple random sample
blvr_df <- data.frame(blvr)
covlist <- c("age", "sex", "race", "crf", "insurance", "income", "hospital_size", "teaching", "region", "ownership", "year", "ps", "linps")
bal.wts1 <- dx.wts(x=blvr_df$wts1, data=blvr_df, vars=covlist,
treat.var="vol_reduction", estimand="ATT")
bal.wts1 type n.treat n.ctrl ess.treat ess.ctrl max.es mean.es max.ks
1 unw 253 265 253 265.0000 2.0306797 0.4589984 0.6194496
2 253 265 253 54.7126 0.3077257 0.1118164 0.2134387
mean.ks iter
1 0.14873095 NA
2 0.05053566 NAThe standardized differences were compared before and after weighting using a Love.plot. The weighting process showed a more satisfactory balance of standardized differences in comparison to the unweighted standardized differences.
bal.table(bal.wts1)$unw
tx.mn tx.sd ct.mn ct.sd std.eff.sz stat
age 67.881 8.390 56.400 15.065 1.368 10.789
sex:Male 0.565 0.496 0.713 0.452 -0.299 12.293
sex:Female 0.435 0.496 0.287 0.452 0.299 NA
race:White 0.901 0.298 0.770 0.421 0.440 8.925
race:Black 0.067 0.250 0.117 0.321 -0.199 NA
race:Other 0.032 0.175 0.113 0.317 -0.466 NA
crf:No 0.538 0.499 0.872 0.334 -0.670 69.838
crf:Yes 0.462 0.499 0.128 0.334 0.670 NA
insurance:Medicare 0.775 0.418 0.408 0.491 0.879 37.291
insurance:Medicaid 0.043 0.204 0.192 0.394 -0.731 NA
insurance:Private 0.182 0.386 0.400 0.490 -0.566 NA
income:0-25th percentile 0.269 0.443 0.279 0.449 -0.024 0.225
income:26th to 50th percentile 0.269 0.443 0.257 0.437 0.027 NA
income:51st to 75th percentile 0.241 0.428 0.264 0.441 -0.054 NA
income:76th to 100th percentile 0.221 0.415 0.200 0.400 0.051 NA
hospital_size:Small 0.107 0.309 0.045 0.208 0.199 3.868
hospital_size:Medium 0.166 0.372 0.208 0.406 -0.112 NA
hospital_size:Large 0.727 0.445 0.747 0.435 -0.045 NA
teaching:Rural 0.028 0.164 0.026 0.160 0.008 0.965
teaching:Urban nonteaching 0.111 0.314 0.075 0.264 0.112 NA
teaching:Urban teaching 0.862 0.345 0.898 0.302 -0.106 NA
region:Northeast 0.126 0.332 0.238 0.426 -0.335 4.752
region:Midwest 0.304 0.460 0.298 0.457 0.014 NA
region:South 0.415 0.493 0.294 0.456 0.245 NA
region:West 0.154 0.361 0.170 0.375 -0.043 NA
ownership:Government 0.123 0.328 0.140 0.347 -0.052 2.065
ownership:Private, non-profit 0.846 0.361 0.792 0.406 0.148 NA
ownership:Private, for-profit 0.032 0.175 0.068 0.252 -0.207 NA
year:2018 0.024 0.152 0.321 0.467 -1.952 50.600
year:2019 0.300 0.458 0.366 0.482 -0.143 NA
year:2020 0.676 0.468 0.313 0.464 0.775 NA
ps 0.719 0.240 0.268 0.258 1.879 20.616
linps 1.288 1.526 -1.811 2.145 2.031 19.031
p ks ks.pval
age 0.000 0.359 0.000
sex:Male 0.000 0.148 0.000
sex:Female NA 0.148 0.000
race:White 0.000 0.131 0.000
race:Black NA 0.050 0.000
race:Other NA 0.082 0.000
crf:No 0.000 0.334 0.000
crf:Yes NA 0.334 0.000
insurance:Medicare 0.000 0.367 0.000
insurance:Medicaid NA 0.149 0.000
insurance:Private NA 0.218 0.000
income:0-25th percentile 0.879 0.010 0.879
income:26th to 50th percentile NA 0.012 0.879
income:51st to 75th percentile NA 0.023 0.879
income:76th to 100th percentile NA 0.021 0.879
hospital_size:Small 0.021 0.061 0.021
hospital_size:Medium NA 0.042 0.021
hospital_size:Large NA 0.020 0.021
teaching:Rural 0.381 0.001 0.381
teaching:Urban nonteaching NA 0.035 0.381
teaching:Urban teaching NA 0.036 0.381
region:Northeast 0.003 0.111 0.003
region:Midwest NA 0.006 0.003
region:South NA 0.121 0.003
region:West NA 0.016 0.003
ownership:Government 0.127 0.017 0.127
ownership:Private, non-profit NA 0.053 0.127
ownership:Private, for-profit NA 0.036 0.127
year:2018 0.000 0.297 0.000
year:2019 NA 0.066 0.000
year:2020 NA 0.363 0.000
ps 0.000 0.619 0.000
linps 0.000 0.619 0.000
[[2]]
tx.mn tx.sd ct.mn ct.sd std.eff.sz stat
age 67.881 8.390 66.767 8.471 0.133 0.923
sex:Male 0.565 0.496 0.613 0.487 -0.096 0.372
sex:Female 0.435 0.496 0.387 0.487 0.096 NA
race:White 0.901 0.298 0.866 0.341 0.118 0.576
race:Black 0.067 0.250 0.104 0.306 -0.148 NA
race:Other 0.032 0.175 0.030 0.170 0.010 NA
crf:No 0.538 0.499 0.655 0.475 -0.236 2.117
crf:Yes 0.462 0.499 0.345 0.475 0.236 NA
insurance:Medicare 0.775 0.418 0.782 0.413 -0.018 0.111
insurance:Medicaid 0.043 0.204 0.034 0.182 0.046 NA
insurance:Private 0.182 0.386 0.184 0.387 -0.005 NA
income:0-25th percentile 0.269 0.443 0.325 0.469 -0.128 0.790
income:26th to 50th percentile 0.269 0.443 0.180 0.384 0.200 NA
income:51st to 75th percentile 0.241 0.428 0.265 0.442 -0.057 NA
income:76th to 100th percentile 0.221 0.415 0.229 0.420 -0.019 NA
hospital_size:Small 0.107 0.309 0.101 0.301 0.019 0.006
hospital_size:Medium 0.166 0.372 0.169 0.375 -0.009 NA
hospital_size:Large 0.727 0.445 0.730 0.444 -0.006 NA
teaching:Rural 0.028 0.164 0.017 0.130 0.063 3.263
teaching:Urban nonteaching 0.111 0.314 0.037 0.188 0.236 NA
teaching:Urban teaching 0.862 0.345 0.946 0.226 -0.245 NA
region:Northeast 0.126 0.332 0.208 0.406 -0.247 1.450
region:Midwest 0.304 0.460 0.334 0.472 -0.065 NA
region:South 0.415 0.493 0.363 0.481 0.105 NA
region:West 0.154 0.361 0.094 0.292 0.166 NA
ownership:Government 0.123 0.328 0.099 0.299 0.071 1.134
ownership:Private, non-profit 0.846 0.361 0.887 0.316 -0.115 NA
ownership:Private, for-profit 0.032 0.175 0.013 0.115 0.104 NA
year:2018 0.024 0.152 0.032 0.176 -0.055 0.187
year:2019 0.300 0.458 0.318 0.466 -0.039 NA
year:2020 0.676 0.468 0.650 0.477 0.056 NA
ps 0.719 0.240 0.663 0.223 0.233 1.767
linps 1.288 1.526 0.818 1.184 0.308 2.489
p ks ks.pval
age 0.356 0.099 0.734
sex:Male 0.542 0.048 0.542
sex:Female NA 0.048 0.542
race:White 0.552 0.035 0.552
race:Black NA 0.037 0.552
race:Other NA 0.002 0.552
crf:No 0.146 0.118 0.146
crf:Yes NA 0.118 0.146
insurance:Medicare 0.870 0.008 0.870
insurance:Medicaid NA 0.009 0.870
insurance:Private NA 0.002 0.870
income:0-25th percentile 0.494 0.057 0.494
income:26th to 50th percentile NA 0.089 0.494
income:51st to 75th percentile NA 0.024 0.494
income:76th to 100th percentile NA 0.008 0.494
hospital_size:Small 0.990 0.006 0.990
hospital_size:Medium NA 0.003 0.990
hospital_size:Large NA 0.003 0.990
teaching:Rural 0.044 0.010 0.044
teaching:Urban nonteaching NA 0.074 0.044
teaching:Urban teaching NA 0.084 0.044
region:Northeast 0.229 0.082 0.229
region:Midwest NA 0.030 0.229
region:South NA 0.052 0.229
region:West NA 0.060 0.229
ownership:Government 0.319 0.023 0.319
ownership:Private, non-profit NA 0.042 0.319
ownership:Private, for-profit NA 0.018 0.319
year:2018 0.784 0.008 0.784
year:2019 NA 0.018 0.784
year:2020 NA 0.026 0.784
ps 0.078 0.213 0.029
linps 0.013 0.213 0.029bal.before.wts1 <- bal.table(bal.wts1)[1]
bal.after.wts1 <- bal.table(bal.wts1)[2]
balance.att.weights <- data.frame(names = rownames(bal.before.wts1$unw),
pre.weighting = 100*bal.before.wts1$unw$std.eff.sz,
ATT.weighted = 100*bal.after.wts1[[1]]$std.eff.sz)
balance.att.weights <- gather(balance.att.weights, timing, szd, 2:3)
balance.att.weights names timing szd
1 age pre.weighting 136.8
2 sex:Male pre.weighting -29.9
3 sex:Female pre.weighting 29.9
4 race:White pre.weighting 44.0
5 race:Black pre.weighting -19.9
6 race:Other pre.weighting -46.6
7 crf:No pre.weighting -67.0
8 crf:Yes pre.weighting 67.0
9 insurance:Medicare pre.weighting 87.9
10 insurance:Medicaid pre.weighting -73.1
11 insurance:Private pre.weighting -56.6
12 income:0-25th percentile pre.weighting -2.4
13 income:26th to 50th percentile pre.weighting 2.7
14 income:51st to 75th percentile pre.weighting -5.4
15 income:76th to 100th percentile pre.weighting 5.1
16 hospital_size:Small pre.weighting 19.9
17 hospital_size:Medium pre.weighting -11.2
18 hospital_size:Large pre.weighting -4.5
19 teaching:Rural pre.weighting 0.8
20 teaching:Urban nonteaching pre.weighting 11.2
21 teaching:Urban teaching pre.weighting -10.6
22 region:Northeast pre.weighting -33.5
23 region:Midwest pre.weighting 1.4
24 region:South pre.weighting 24.5
25 region:West pre.weighting -4.3
26 ownership:Government pre.weighting -5.2
27 ownership:Private, non-profit pre.weighting 14.8
28 ownership:Private, for-profit pre.weighting -20.7
29 year:2018 pre.weighting -195.2
30 year:2019 pre.weighting -14.3
31 year:2020 pre.weighting 77.5
32 ps pre.weighting 187.9
33 linps pre.weighting 203.1
34 age ATT.weighted 13.3
35 sex:Male ATT.weighted -9.6
36 sex:Female ATT.weighted 9.6
37 race:White ATT.weighted 11.8
38 race:Black ATT.weighted -14.8
39 race:Other ATT.weighted 1.0
40 crf:No ATT.weighted -23.6
41 crf:Yes ATT.weighted 23.6
42 insurance:Medicare ATT.weighted -1.8
43 insurance:Medicaid ATT.weighted 4.6
44 insurance:Private ATT.weighted -0.5
45 income:0-25th percentile ATT.weighted -12.8
46 income:26th to 50th percentile ATT.weighted 20.0
47 income:51st to 75th percentile ATT.weighted -5.7
48 income:76th to 100th percentile ATT.weighted -1.9
49 hospital_size:Small ATT.weighted 1.9
50 hospital_size:Medium ATT.weighted -0.9
51 hospital_size:Large ATT.weighted -0.6
52 teaching:Rural ATT.weighted 6.3
53 teaching:Urban nonteaching ATT.weighted 23.6
54 teaching:Urban teaching ATT.weighted -24.5
55 region:Northeast ATT.weighted -24.7
56 region:Midwest ATT.weighted -6.5
57 region:South ATT.weighted 10.5
58 region:West ATT.weighted 16.6
59 ownership:Government ATT.weighted 7.1
60 ownership:Private, non-profit ATT.weighted -11.5
61 ownership:Private, for-profit ATT.weighted 10.4
62 year:2018 ATT.weighted -5.5
63 year:2019 ATT.weighted -3.9
64 year:2020 ATT.weighted 5.6
65 ps ATT.weighted 23.3
66 linps ATT.weighted 30.8ggplot(balance.att.weights, aes(x = szd, y = reorder(names, szd), color = timing)) +
geom_point(size = 2, alpha = 0.5) +
geom_vline(xintercept = 0) +
geom_vline(xintercept = c(-10,10), linetype = "dashed", col = "blue") +
labs(x = "Standardized Differences", y = "",
title = "Love Plot of Standarized Diffs in The BLVR Data",
subtitle = "Before and After ATT Weighting") 
Rubin’s rules were obtained after weighting. Rubin’s rule 1 was obtained by using the standardized effect size after weighting for the linear propensity score (0.28). When multiplied by 100, the result was 28% indicating that Rubin’s rule 1 was satisfied (although not perfectly). Rubin’s rule 2 was calculated by dividing the square root of treatment group standard deviation by the square root of control group standard deviation, resulting in a value of 1.54. Rubin’s rule 2 was also satisfied, although less perfectly.
arfblvr.design <- svydesign(ids=~1, weights=~wts1, data=blvr)weighted_arf <- svyglm(arf ~ vol_reduction, design=blvr.design, family=quasibinomial())kable(tidy(weighted_arf, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_material() | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.13 | 0.32 | -6.49 | 0.00 | 0.07 | 0.24 |
| vol_reduction | 1.68 | 0.36 | 1.43 | 0.15 | 0.83 | 3.40 |
The estimate of treatment effect on acute respiratory failure was obtained by a regression model using survey weighting and a quasi binomial design.
The weighted odds ratio estimate for the BLVR group was 1.68, with a 95% confidence interval (0.83,3.4). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.68 the odds of respiratory failure compared to the person receiving LVRS. It should be noted that this confidence interveal crosses 1
ptxweighted_ptx <- svyglm(ptx ~ vol_reduction, design=blvr.design, family=quasibinomial())kable(tidy(weighted_ptx, conf.int = TRUE, conf.level = 0.95, exponentiate = TRUE), digits = 2) |> kable_material() | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.73 | 0.28 | -1.13 | 0.26 | 0.43 | 1.26 |
| vol_reduction | 0.36 | 0.32 | -3.21 | 0.00 | 0.19 | 0.67 |
The weighted odds ratio estimate for the `BLVR group was 0.36, with a 95% confidence interval (0.19,0.67). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.36 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1
A logistic regression for the outcomes of acute respiratory failure and pneumothorax was fitted with a regression adjustment for the linear propensity score to obtain a doubly-robust estimate and account for residual imbalance after weighting
robust.design <- svydesign(ids=~1, weights=~wts1, data=blvr)robust_arf <- svyglm(arf ~ vol_reduction + linps, design=robust.design, family=quasibinomial())kable(tidy(robust_arf, exponentiate = TRUE, conf.int = TRUE, conf.level = 0.95), digits = 2) |> kable_material()| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.15 | 0.30 | -6.23 | 0.00 | 0.08 | 0.28 |
| vol_reduction | 1.89 | 0.36 | 1.76 | 0.08 | 0.93 | 3.87 |
| linps | 0.75 | 0.08 | -3.58 | 0.00 | 0.64 | 0.88 |
The doubly-robust odds ratio estimate for the BLVR group was 1.89, with a 95% confidence interval (0.93,3.87). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.89 the odds of acute respiratory failure compared to the person receiving LVRS. The confidence interval here crosses 1
ptxrobust_ptx <- svyglm(ptx ~ vol_reduction + linps, design=robust.design, family=quasibinomial())And then the Odds Ratio with the 95% confidence interval
kable(tidy(robust_ptx, conf.int = TRUE, exponentiate = TRUE, conf.level = 0.95), digits = 2) |> kable_material()| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | 0.83 | 0.23 | -0.83 | 0.41 | 0.53 | 1.29 |
| vol_reduction | 0.38 | 0.34 | -2.79 | 0.01 | 0.20 | 0.75 |
| linps | 0.86 | 0.12 | -1.30 | 0.20 | 0.68 | 1.08 |
The doubly-robust odds ratio estimate for the `BLVR group was 0.38, with a 95% confidence interval (0.2,0.75). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.38 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1
arfarf_blvr <- tibble(
analysis = c("Unadjusted", "Matched",
"ATT Weighted", "Double-Robust"),
estimate = c(1.65, 1.73, 1.68, 1.89),
conf.low = c(1, 0.82, 0.83, 0.93),
conf.high = c(2.74, 3.63, 3.4, 3.87))
ggplot(arf_blvr, aes(x = analysis, y = estimate)) +
geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = 0.5, color = "orange3") +
geom_label(aes(label = estimate), size = 6, color = "white", fill = "lightblue3") +
labs(title = "Comparison Of Odds Ratio Estimates", subtitle = "Between Different Adjustment Approaches In Respiratory Failure", x = "Method of PS Adjustment", y = "Odds Ratio")
| Method | Estimate | Lower 95% CI | upper 95% CI |
|---|---|---|---|
| Unadjusted | 1.65 | 1.00 | 2.74 |
| Matched | 1.73 | 0.82 | 3.63 |
| Weighted | 1.68 | 0.83 | 3.4 |
| Double-Robust | 1.89 | 0.93 | 3.87 |
Based on the comparison above, the odds ratio of acute respiratory failure didn’t substantially change from its estimate prior to propensity score adjustment. All of the confidence intervals contain 1.
ptxptx_blvr <- tibble(
analysis = c("Unadjusted", "Matched",
"ATT Weighted", "Double-Robust"),
estimate = c(0.4, 0.68, 0.36, 0.38),
conf.low = c(0.27, 0.38, 0.19, 0.2),
conf.high = c(0.58, 1.22, 0.67, 0.75))
ggplot(ptx_blvr, aes(x = analysis, y = estimate)) +
geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width = 0.5, color = "purple3") +
geom_label(aes(label = estimate), size = 6, color = "white", fill = "lightpink3") +
labs(title = "Comparison Of Odds Ratio Estimates", subtitle = "Between Different Adjustment Approaches In Pneumothorax", x = "Method of PS Adjustment", y = "Odds Ratio")
| Method | Estimate | Lower 95% CI | upper 95% CI |
|---|---|---|---|
| Unadjusted | 0.4 | 0.27 | 0.58 |
| Matched | 0.68 | 0.38 | 1.22 |
| Weighted | 0.36 | 0.19 | 0.67 |
| Double-Robust | 0.38 | 0.2 | 0.75 |
The comparison above showed that the odds ratio for pneumothorax also didn’t substantially change after propensity score adjustment, reflecting a plausible conclusion that patients undergoing lung volume reduction surgery tended to have higher odds of pneumothorax based on the nature of invasive surgical intervention on the lung.
Is bronchoscopic lung volume reduction associated with higher odds of developing post-procedural acute respiratory failure in adult patients with emphysema compared to lung volume reduction surgery?
Yes. Based on the odds ratio obtained prior to propensity score adjustment and after, the patients undergoing bronchoscopic lung volume reduction surgery tended to have higher odds of developing acute respiratory failure, granted that the confidence interval for this estimate always crossed 1
It was noted that despite propensity score matching, weighting and a double robust adjustment (Section 19.1), the odds ratio of acute respiratory failure did not show a meaningful change from prior to adjustment. The small fluctuations wouldn’t have significant clinical implication.
Even though the confidence intervals crossed 1, this analysis showed an interesting trend in the odds of respiratory failure. A possible explanation could be that there was a higher percentage of patients with chronic respiratory failure undergoing bronchoscopic lung volume reduction in comparison to patients undergoing surgical lung volume reduction, so the patient’s lungs have already been compromised and further interventions
Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural pneumothorax in adult patients with emphysema compared to lung volume reduction surgery?
No. BLVR was associated with lower odds of pneumothorax (Section 19.2) despite propensity score adjustment methods, which was rather plausible, as the surgical procedure was more invasive and expected to cause a higher incidence of pneumothorax. The presence of pneumothorax is well described in the literature, and could indicate a favorable response after lung volume reduction
ggplot(blvr, aes(x = year, fill = vol_reduction_f)) +
geom_bar(position = "dodge", col = "white") +
scale_fill_viridis_d(option = "turbo", alpha = 0.5) +
labs(title = "Lung Volume Reduction", subtitle = "Distribution By Year", x = "Volume Reduction", y = "# of Patients", fill = "Volume Reduction") 
Bronchoscopic lung volume reduction with endobronchial valves has largely replaced lung volume reduction surgery as a possible solution to improve excessive dyspnea. Patient selection is a key element in developing post procedural complications that could be life threatening.
Understanding the potential complications that may arise in real-world settings outside controlled clinical trials becomes increasingly crucial as the number of procedures grows and their availability expands beyond large referral centers. This knowledge aids physicians in better patient selection for the procedure, enabling them to implement effective action plans. For instance, the introduction of a “Pneumothorax Travel Cart” containing procedural kits for pneumothorax management has significantly improved time to management of pneumothorax after bronchoscopic lung volume reduction.
By trying to get a sense of the complication rates nationally, this information could be used to obtain further data from endobronchial valve device companies. These data could include more patient covariates that play a role in procedural selection, and provide a better characterization of the complications encountered.
By doing this project, I got a chance to apply all the matching methods (Section 14) we learned over the span of this course. While I was slightly distraught at the beginning when I ended up with a sample with more treated patients than controls, I used it as an opportunity to go over all the plausible methods to properly troubleshoot and be prepared for similar scenarios in the future beyond the security of this class.
xfun::session_info()R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22631)
Locale:
LC_COLLATE=English_United States.utf8
LC_CTYPE=English_United States.utf8
LC_MONETARY=English_United States.utf8
LC_NUMERIC=C
LC_TIME=English_United States.utf8
Package version:
abind_1.4-5 askpass_1.2.0 backports_1.4.1
base64enc_0.1-3 bit_4.0.5 bit64_4.0.5
blob_1.2.4 boot_1.3-28.1 brio_1.1.4
broom_1.0.5 bslib_0.6.1 cachem_1.0.8
callr_3.7.3 car_3.1-2 carData_3.0-5
cellranger_1.1.0 chk_0.9.1 class_7.3-22
cli_3.6.2 clipr_0.8.0 cobalt_4.5.3
codetools_0.2-19 colorspace_2.1-0 compiler_4.3.2
conflicted_1.2.0 cpp11_0.4.7 crayon_1.5.2
curl_5.2.0 data.table_1.14.10 DBI_1.2.1
dbplyr_2.4.0 deldir_2.0-2 desc_1.4.3
diffobj_0.3.5 digest_0.6.34 dplyr_1.1.4
dtplyr_1.3.1 e1071_1.7-14 ellipsis_0.3.2
evaluate_0.23 fansi_1.0.6 farver_2.1.1
fastmap_1.1.1 fontawesome_0.5.2 forcats_1.0.0
foreach_1.5.2 fs_1.6.3 gargle_1.5.2
gbm_2.1.9 gdata_3.0.0 generics_0.1.3
ggformula_0.12.0 ggplot2_3.4.4 ggridges_0.5.5
glmnet_4.1-8 glue_1.7.0 gmodels_2.18.1.1
googledrive_2.1.1 googlesheets4_1.1.1 graphics_4.3.2
grDevices_4.3.2 grid_4.3.2 gridExtra_2.3
gtable_0.3.4 gtools_3.9.5 haven_2.5.4
highr_0.10 hms_1.1.3 htmltools_0.5.7
htmlwidgets_1.6.4 httr_1.4.7 ids_1.0.1
interp_1.1-6 isoband_0.2.7 iterators_1.0.14
janitor_2.2.0 jomo_2.7-6 jpeg_0.1-10
jquerylib_0.1.4 jsonlite_1.8.8 kableExtra_1.3.4
knitr_1.45 labeling_0.4.3 labelled_2.12.0
lattice_0.22-5 latticeExtra_0.6-30 lifecycle_1.0.4
lme4_1.1-35.1 lubridate_1.9.3 magrittr_2.0.3
MASS_7.3-60.0.1 Matching_4.10-14 Matrix_1.6-1.1
MatrixModels_0.5-3 memoise_2.0.1 methods_4.3.2
mgcv_1.9.0 mice_3.16.0 mime_0.12
minqa_1.2.6 mitml_0.4-5 mitools_2.4
modelr_0.1.11 mosaic_1.9.0 mosaicCore_0.9.4.0
mosaicData_0.20.4 munsell_0.5.0 naniar_1.0.0
nlme_3.1-163 nloptr_2.0.3 nnet_7.3-19
norm_1.0.11.1 numDeriv_2016.8.1.1 openssl_2.1.1
ordinal_2023.12.4 pan_1.9 parallel_4.3.2
pbkrtest_0.5.2 pillar_1.9.0 pkgbuild_1.4.3
pkgconfig_2.0.3 pkgload_1.3.3 plyr_1.8.9
png_0.1-8 praise_1.0.0 prettyunits_1.2.0
processx_3.8.3 progress_1.2.3 proxy_0.4-27
ps_1.7.5 purrr_1.0.2 quantreg_5.97
R6_2.5.1 ragg_1.2.7 rappdirs_0.3.3
RColorBrewer_1.1-3 Rcpp_1.0.12 RcppEigen_0.3.4.0.0
readr_2.1.5 readxl_1.4.3 rematch_2.0.0
rematch2_2.1.2 repr_1.1.6 reprex_2.1.0
rlang_1.1.3 rmarkdown_2.25 rpart_4.1.23
rprojroot_2.0.4 rstudioapi_0.15.0 rvest_1.0.3
sass_0.4.8 scales_1.3.0 selectr_0.4.2
shape_1.4.6 skimr_2.1.5 snakecase_0.11.1
SparseM_1.81 splines_4.3.2 stats_4.3.2
stringi_1.8.3 stringr_1.5.1 survey_4.2-1
survival_3.5-7 svglite_2.1.3 sys_3.4.2
systemfonts_1.0.5 tableone_0.13.2 testthat_3.2.1
textshaping_0.3.7 tibble_3.2.1 tidyr_1.3.0
tidyselect_1.2.0 tidyverse_2.0.0 timechange_0.2.0
tinytex_0.49 tools_4.3.2 twang_2.6
tzdb_0.4.0 ucminf_1.2.1 UpSetR_1.4.0
utf8_1.2.4 utils_4.3.2 uuid_1.1.1
vctrs_0.6.5 viridis_0.6.4 viridisLite_0.4.2
visdat_0.6.0 vroom_1.6.5 waldo_0.5.2
webshot_0.5.5 withr_2.5.2 xfun_0.41
xgboost_1.7.7.1 xml2_1.3.6 xtable_1.8-4
yaml_2.3.8 zoo_1.8-12